library(dplyr) # ungroup()
library(ggtree) # BuildClusterTree()
library(gridExtra) # grid.arrange()
library(gtools) # smartbind()
library(parallel) # detectCores()
library(plotly) # plot_ly()
library(Seurat) # Idents()
library(SeuratWrappers) # RunPrestoAll()
library(ShinyCell) # createConfig()
library(tidyr) # %>%
out <- "../../results/all_clusters/"
sample_order <- c("E3.M","E3.F","E4.M","E4.F")
sample_colors <- c("#26946A","#1814A1","#EAC941","#DF5F00")
sample_order2 <- c("Male E3", "Male E4", "Female E3", "Female E4")
isoform_order <- c("E4","E3")
isoform_colors <- c("darkgray","cornflowerblue")
sex_order <- c("Male","Female")
sex_colors <- c("darkgray","purple")
mouse.unannotated <- readRDS("../../rObjects/unannotated_obj.rds")
Idents(mouse.unannotated) <- "seurat_clusters"
DefaultAssay(mouse.unannotated) <- "RNA"
mouse.unannotated <- NormalizeData(mouse.unannotated)
## Normalizing layer: counts
mouse.unannotated
## An object of class Seurat
## 41866 features across 16965 samples within 2 assays
## Active assay: RNA (21089 features, 0 variable features)
## 2 layers present: data, counts
## 1 other assay present: SCT
## 3 dimensional reductions calculated: pca, harmony, umap
cluster_colors <- c("red4","red","orange","gold","yellow","lightgreen","chartreuse3","darkgreen",
"cyan","steelblue","blue","deeppink","pink","salmon","purple3","orchid","tan","chocolate4",
"gray30","gray70","black","aquamarine")
u1 <- DimPlot(object = mouse.unannotated,
reduction = "umap",
shuffle = TRUE,
repel = TRUE,
raster = FALSE,
cols = cluster_colors,
label = TRUE)
u1
u2 <- DimPlot(object = mouse.unannotated,
reduction = "umap",
shuffle = TRUE,
repel = TRUE,
dims = c(2,3),
cols = cluster_colors,
label = TRUE)
u2
# UMAP percent.mt
FeaturePlot(mouse.unannotated,
reduction = "umap",
features = "percent.mt") +
scale_colour_gradientn(colours = c("blue","lightblue","yellow","orange","red"))
# UMAP percent.ribo
FeaturePlot(mouse.unannotated,
reduction = "umap",
features = "percent.ribo.protein") +
scale_colour_gradientn(colours = c("blue","lightblue","yellow","orange","red"))
# UMAP percent.hb
FeaturePlot(mouse.unannotated,
reduction = "umap",
features = "percent.hb") +
scale_colour_gradientn(colours = c("blue","lightblue","yellow","orange","red"))
# UMAP nCount
FeaturePlot(mouse.unannotated,
reduction = "umap",
features = "nCount_RNA") +
scale_colour_gradientn(colours = c("blue","lightblue","yellow","orange","red"))
# UMAP nFeature
FeaturePlot(mouse.unannotated,
reduction = "umap",
features = "nFeature_RNA") +
scale_colour_gradientn(colours = c("blue","lightblue","yellow","orange","red"))
# UMAP cell.complexity
FeaturePlot(mouse.unannotated,
reduction = "umap",
features = "cell.complexity") +
scale_colour_gradientn(colours = c("blue","lightblue","yellow","orange","red"))
VlnPlot(mouse.unannotated,
features = "nCount_RNA",
split.by = "seurat_clusters")
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##
## This message will be shown once per session.
VlnPlot(mouse.unannotated,
features = "nFeature_RNA",
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "cell.complexity",
split.by = "seurat_clusters")
# Cells per sample per cluster
sample_ncells <- FetchData(mouse.unannotated,
vars = c("ident", "sample")) %>%
dplyr::count(ident,sample) %>%
tidyr::spread(ident, n)
sample_ncells
## sample 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## 1 E3.M 598 529 228 220 267 254 189 282 344 196 118 110 84 105 94 49 29 94
## 2 E3.F 288 578 302 326 199 229 184 285 27 116 178 60 65 31 40 21 40 58
## 3 E4.M 581 381 238 218 232 242 228 304 313 238 165 113 63 66 74 55 37 48
## 4 E4.F 818 769 596 425 366 312 419 109 289 337 295 151 155 142 122 182 171 70
## 18 19 20 21
## 1 41 64 59 18
## 2 70 7 43 8
## 3 62 80 62 21
## 4 82 93 59 55
# Cells per isoform per cluster
isoform_ncells <- FetchData(mouse.unannotated,
vars = c("ident", "isoform")) %>%
dplyr::count(ident,isoform) %>%
tidyr::spread(ident, n)
isoform_ncells
## isoform 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 1 E4 1399 1150 834 643 598 554 647 413 602 575 460 264 218 208 196 237 208
## 2 E3 886 1107 530 546 466 483 373 567 371 312 296 170 149 136 134 70 69
## 17 18 19 20 21
## 1 118 144 173 121 76
## 2 152 111 71 102 26
# Cells per sex per cluster
sex_ncells <- FetchData(mouse.unannotated,
vars = c("ident", "sex")) %>%
dplyr::count(ident,sex) %>%
tidyr::spread(ident, n)
sex_ncells
## sex 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
## 1 Male 1179 910 466 438 499 496 417 586 657 434 283 223 147 171 168 104 66
## 2 Female 1106 1347 898 751 565 541 603 394 316 453 473 211 220 173 162 203 211
## 17 18 19 20 21
## 1 142 103 144 121 39
## 2 128 152 100 102 63
# User params
goi <- "Malat1"
threshold <- 1
# Subset data
log2.threshold <- log2(threshold + 0.01)
counts.df <- FetchData(mouse.unannotated, vars = goi)
colnames(counts.df) <- "counts"
log2.counts.df <- log2(counts.df + 0.01)
# Histogram
title <- paste0(goi, "\nnCount_RNA > ", threshold)
hist1 <- ggplot(counts.df, aes(x = counts)) +
geom_histogram(bins = 100, fill = "gray", color = "black") +
labs(title = title, x=NULL, y=NULL) +
xlab(paste0(goi, " nCount_RNA")) + ylab("# of Samples") + theme_bw() +
geom_vline(xintercept = threshold, col = "blue") +
annotate("rect",
xmin = -Inf,
xmax = threshold,
ymin = 0,
ymax=Inf,
alpha=0.2,
fill="chocolate4") +
annotate("rect",
xmin = threshold,
xmax = Inf,
ymin = 0,
ymax=Inf,
alpha=0.2,
fill="deepskyblue")
# Histogram log transformed
hist2 <- ggplot(log2.counts.df, aes(x = counts)) +
geom_histogram(bins = 100, fill = "gray", color = "black") +
labs(title = title, x=NULL, y=NULL) +
xlab(paste0("Log2(",goi, " nCount_RNA)")) + ylab("# of Samples") + theme_bw() +
geom_vline(xintercept = log2.threshold, col = "blue") +
annotate("rect",
xmin = -Inf,
xmax = log2.threshold,
ymin = 0,
ymax=Inf,
alpha=0.2,
fill="chocolate4") +
annotate("rect",
xmin = log2.threshold,
xmax = Inf,
ymin = 0,
ymax=Inf,
alpha=0.2,
fill="deepskyblue")
# plot
plots1 <- list(hist1,hist2)
layout1 <- rbind(c(1),c(2))
grid1 <- grid.arrange(grobs = plots1, layout_matrix = layout1)
# user define variable
goi <- "Malat1"
# Extract counts data
DefaultAssay(mouse.unannotated) <- "RNA"
Idents(mouse.unannotated) <- "seurat_clusters"
geneData <- FetchData(mouse.unannotated,
vars = goi,
slot = "counts")
## Warning: The `slot` argument of `FetchData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
geneData <- geneData > 1
table(geneData)
## geneData
## FALSE TRUE
## 1248 15717
mouse.unannotated$Expression <- geneData
FetchData(mouse.unannotated,
vars = c("ident", "Expression")) %>%
dplyr::count(ident, Expression) %>%
tidyr::spread(ident, n)
## Expression 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
## 1 FALSE 2 3 4 2 43 1 690 9 NA NA NA NA 285 13 122
## 2 TRUE 2283 2254 1360 1187 1021 1036 330 971 973 887 756 434 82 331 208
## 15 16 17 18 19 20 21
## 1 2 65 1 4 1 1 NA
## 2 305 212 269 251 243 222 102
# Plot
mouse.unannotated@meta.data %>%
group_by(seurat_clusters, Expression) %>%
dplyr::count() %>%
group_by(seurat_clusters) %>%
dplyr::mutate(percent = 100*n/sum(n)) %>%
ungroup() %>%
ggplot(aes(x=seurat_clusters,y=percent, fill=Expression)) +
geom_col() +
ggtitle(paste0("Percentage of cells with > 1 counts for ", goi)) +
theme(axis.text.x = element_text(angle = 90))
mouse.unannotated <- BuildClusterTree(object = mouse.unannotated,
dims = 1:15,
reorder = FALSE,
reorder.numeric = FALSE)
tree <- mouse.unannotated@tools$BuildClusterTree
tree$tip.label <- paste0(tree$tip.label)
ggtree::ggtree(tree, aes(x, y)) +
scale_y_reverse() +
ggtree::geom_tree() +
ggtree::theme_tree() +
ggtree::geom_tiplab(offset = 1) +
ggtree::geom_tippoint(color = cluster_colors[1:length(tree$tip.label)], shape = 16, size = 5) +
coord_cartesian(clip = 'off') +
theme(plot.margin = unit(c(0,2.5,0,0), 'cm'))
VlnPlot(mouse.unannotated,
features = "Ptprc",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Lyve1",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Pecam1",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Cd19",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Fcrla",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Cd79a",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Cd79b",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Sdc1",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Trac",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Cd3e",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Cd8a",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Cd4",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Ly6c1",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Flt1",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Col1a1",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Col1a2",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Fcer1a",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Kit", # aka Cd117
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "C1qa",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "C1qb",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = c("Itgax","Cd209a","Flt3","Zbtb46","Ccr2"),
cols = cluster_colors,
split.by = "seurat_clusters",
stack = TRUE,
flip = TRUE)
VlnPlot(mouse.unannotated,
features = "Ly6g",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Retnlg",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Acta2",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Myl9",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Rgs5",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Cdh19",
cols = cluster_colors,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = "Mpz",
cols = cluster_colors,
split.by = "seurat_clusters")
Idents(mouse.unannotated) <- "seurat_clusters"
goi <- c("Pdgfrb","Acta2","Vwf","Cldn5","Pecam1","Flt4","Prox1","Lyve1","Ptprc",
"Cd19","Ms4a1","Ighd","Igha","Sdc1","Cd3e","Trbc2","Il7r","Nkg7","Klrb1b","Klrb1c",
"Gata3","Rora","Itgax","H2-Eb1","Ccr2","Ly6c2","Lyz2","Ly6g","Itgam",
"Mrc1","Csf1r","Cd38","Mki67","Mcpt4","Ms4a2","Col1a2","Plp1")
v1 <- VlnPlot(mouse.unannotated,
features = goi,
cols = cluster_colors,
split.by = "seurat_clusters",
flip = TRUE,
stack = TRUE)
v1
pdf(paste0(out, "markers/sandro_markers_violin_unannotated.pdf"), width = 8, height = 8)
v1
dev.off()
## png
## 2
# auto find markers
Idents(mouse.unannotated) <- "seurat_clusters"
all.markers <- SeuratWrappers::RunPrestoAll(
object = mouse.unannotated,
assay = "RNA",
slot = "counts",
only.pos = TRUE
)
## Calculating cluster 0
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 5
## Calculating cluster 6
## Calculating cluster 7
## Calculating cluster 8
## Calculating cluster 9
## Calculating cluster 10
## Calculating cluster 11
## Calculating cluster 12
## Calculating cluster 13
## Calculating cluster 14
## Calculating cluster 15
## Calculating cluster 16
## Calculating cluster 17
## Calculating cluster 18
## Calculating cluster 19
## Calculating cluster 20
## Calculating cluster 21
# filter based on p_val_adj
all.markers <- all.markers[all.markers$p_val_adj < 0.01,]
# subset table to top 2 and top 20 markers
top2 <- Reduce(rbind,
by(all.markers,
all.markers["cluster"],
head,
n = 2))
top20 <- Reduce(rbind,
by(all.markers,
all.markers["cluster"],
head,
n = 20))
# plot violin
v1 <- VlnPlot(mouse.unannotated,
features = top2$gene,
split.by = "seurat_clusters",
flip = TRUE,
stack = TRUE,
cols = cluster_colors)
v1
# save table
write.table(all.markers,
paste0(out, "markers/unannotated_auto_find_markers_adjpval_0.01.tsv"),
quote = FALSE,
row.names = FALSE)
# compare
table(all.markers$cluster)
##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
## 2243 3010 3430 2888 1618 3451 1322 1304 2491 1382 1344 4527 570 2207 491 867
## 16 17 18 19 20 21
## 536 649 1739 1986 1138 538
# subset based on cluster
cluster0 <- all.markers[all.markers$cluster == 0,]
cluster1 <- all.markers[all.markers$cluster == 1,]
cluster2 <- all.markers[all.markers$cluster == 2,]
cluster3 <- all.markers[all.markers$cluster == 3,]
cluster4 <- all.markers[all.markers$cluster == 4,]
cluster5 <- all.markers[all.markers$cluster == 5,]
cluster6 <- all.markers[all.markers$cluster == 6,]
cluster7 <- all.markers[all.markers$cluster == 7,]
cluster8 <- all.markers[all.markers$cluster == 8,]
cluster9 <- all.markers[all.markers$cluster == 9,]
cluster10 <- all.markers[all.markers$cluster == 10,]
cluster11 <- all.markers[all.markers$cluster == 11,]
cluster12 <- all.markers[all.markers$cluster == 12,]
cluster13 <- all.markers[all.markers$cluster == 13,]
cluster14 <- all.markers[all.markers$cluster == 14,]
cluster15 <- all.markers[all.markers$cluster == 15,]
cluster16 <- all.markers[all.markers$cluster == 16,]
cluster17 <- all.markers[all.markers$cluster == 17,]
cluster18 <- all.markers[all.markers$cluster == 18,]
cluster19 <- all.markers[all.markers$cluster == 19,]
cluster20 <- all.markers[all.markers$cluster == 20,]
cluster21 <- all.markers[all.markers$cluster == 21,]
VlnPlot(mouse.unannotated,
features = cluster0$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster1$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster2$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster3$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster4$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster5$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster6$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster7$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster8$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster9$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster10$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster11$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster12$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster13$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster14$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster15$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster16$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster17$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster18$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster19$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster20$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
VlnPlot(mouse.unannotated,
features = cluster21$gene[1:20],
cols = cluster_colors,
stack = TRUE,
flip = TRUE,
split.by = "seurat_clusters")
# Rename all identities
mouse.annotated <- RenameIdents(object = mouse.unannotated,
"0" = "BECs",
"1" = "Macrophages",
"2" = "Fibroblasts",
"3" = "Macrophages",
"4" = "BECs",
"5" = "Dendritic cells",
"6" = "Macrophages",
"7" = "Neutrophils",
"8" = "B cells",
"9" = "T and NK cells",
"10" = "ILCs",
"11" = "T and NK cells",
"12" = "Fibroblasts",
"13" = "Pericytes and SMCs",
"14" = "BECs",
"15" = "Schwann cells",
"16" = "Fibroblasts",
"17" = "Mast cells",
"18" = "Myeloid precursors",
"19" = "B cells",
"20" = "Plasma cells",
"21" = "LECs")
# save idents
mouse.annotated$annotated_clusters <- Idents(mouse.annotated)
# set levels
mouse.annotated$annotated_clusters <- factor(mouse.annotated$annotated_clusters,
levels = c("Pericytes and SMCs",
"BECs",
"LECs",
"B cells",
"Plasma cells",
"T and NK cells",
"ILCs",
"Dendritic cells",
"Neutrophils",
"Macrophages",
"Myeloid precursors",
"Mast cells",
"Fibroblasts",
"Schwann cells"))
# set ident
Idents(mouse.annotated) <- "annotated_clusters"
# set params
DefaultAssay(mouse.annotated) <- "RNA"
mouse.annotated <- NormalizeData(mouse.annotated)
## Normalizing layer: counts
cluster_colors <- c("#B5B9BA", # Pericytes and SMCs
"#3385BB", # BECs
"#40BBFF", # LECs
"#a5d5a9", # B cells
"#5dbd64", # Plasma cells
"#1c7e24", # T and NK cells
"#F57C7C", # ILCs
"#E42622", # Dendritic cells
"#FBB268", # Neutrophils
"#FE8D19", # Macrophages
"#DE9E83", # Myeloid precursors
"#A6CEE3", # Mast cells
"#9D7BBA", # Fibroblasts
"#977899") # Schwann cells
# umap
umap1 <- DimPlot(object = mouse.annotated,
reduction = "umap",
repel = TRUE,
group.by = "annotated_clusters",
cols = cluster_colors)
umap1
umap2 <- DimPlot(object = mouse.annotated,
reduction = "umap",
repel = TRUE,
dims = c(2,3),
group.by = "annotated_clusters",
cols = cluster_colors)
umap2
# auto find markers
Idents(mouse.annotated) <- "annotated_clusters"
all.markers <- SeuratWrappers::RunPrestoAll(
object = mouse.annotated,
assay = "RNA",
slot = "counts",
only.pos = TRUE
)
## Calculating cluster Pericytes and SMCs
## Calculating cluster BECs
## Calculating cluster LECs
## Calculating cluster B cells
## Calculating cluster Plasma cells
## Calculating cluster T and NK cells
## Calculating cluster ILCs
## Calculating cluster Dendritic cells
## Calculating cluster Neutrophils
## Calculating cluster Macrophages
## Calculating cluster Myeloid precursors
## Calculating cluster Mast cells
## Calculating cluster Fibroblasts
## Calculating cluster Schwann cells
# filter based on p_val_adj
all.markers <- all.markers[all.markers$p_val_adj < 0.01,]
# subset table to top 2 and top 20
top2 <- Reduce(rbind,
by(all.markers,
all.markers["cluster"],
head,
n = 2))
top20 <- Reduce(rbind,
by(all.markers,
all.markers["cluster"],
head,
n = 20))
# save
write.table(all.markers,
paste0(out, "markers/annotated_auto_find_markers_adjpval_0.01.tsv"),
quote = FALSE,
row.names = FALSE)
v1 <- VlnPlot(mouse.annotated,
features = top2$gene,
split.by = "annotated_clusters",
flip = TRUE,
stack = TRUE,
cols = cluster_colors)
v1
goi <- c("Pdgfrb","Acta2","Vwf","Cldn5","Pecam1","Flt4","Prox1","Lyve1","Ptprc",
"Cd19","Ms4a1","Ighd","Igha","Sdc1","Cd3e","Trbc2","Il7r","Nkg7","Klrb1b","Klrb1c",
"Gata3","Rora","Itgax","H2-Eb1","Ccr2","Ly6c2","Lyz2","Ly6g","Itgam",
"Mrc1","Csf1r","Cd38","Mki67","Mcpt4","Ms4a2","Col1a2","Plp1")
v2 <- VlnPlot(mouse.annotated,
group.by = "annotated_clusters",
split.by = "annotated_clusters",
stack = TRUE,
flip = TRUE,
cols = cluster_colors,
features = goi)
v2
# Step 1: Pseudo-bulk the counts based on sample and cell type
pseudo <- AggregateExpression(
object = mouse.annotated,
assays = "RNA",
features = rownames(mouse.annotated),
return.seurat = TRUE,
group.by = c("sample", "annotated_clusters")
)
## Centering and scaling data matrix
# Step 2: Normalize the data
pseudo <- NormalizeData(pseudo,
normalization.method = "LogNormalize",
scale.factor = 10000)
## Normalizing layer: counts
# Step 3: Find variable features
pseudo <- FindVariableFeatures(pseudo,
selection.method = "vst",
nfeatures = 2000)
## Finding variable features for layer counts
# Step 4: Scale the data
pseudo <- ScaleData(pseudo,
features = rownames(pseudo))
## Centering and scaling data matrix
# Step 5: Run PCA
pseudo <- RunPCA(pseudo,
features = VariableFeatures(pseudo),
npcs = 10)
## PC_ 1
## Positive: Dclre1c, Bcl2a1b, Ms4a6b, Fcgr2b, Rgs10, Gpr65, Cd86, Ctsc, Cfp, Gm26542
## Clcn5, Tmem8, Zeb2os, Fcgr3, Gpr160, Slc37a2, Nlrc4, Apobec1, Csf1r, P2ry12
## Cd300c2, Ccr5, Adgre1, Gm21188, Ifi213, Rbpj, Clec10a, Cd68, Cndp2, C1qb
## Negative: Grb10, Plpp3, Fbxl7, Cttn, Ppic, Ghr, Uaca, Smad6, Synpo, Cachd1
## Adarb1, Yap1, Serpinh1, Tmtc1, Bgn, Sulf1, Col5a2, Tcaf1, Wwtr1, Pdgfd
## Fam20a, Prdm5, Prickle1, Nfib, Fat4, Klf12, Mn1, Nr1d1, Id1, Nid1
## PC_ 2
## Positive: Serpinf1, Slc35g1, Olfml3, Clmp, Frmpd4, Shisa3, Car13, H19, Egfr, Col12a1
## Ecrg4, Cdon, Sfrp4, Mfap4, Shank1, Clec11a, Lmx1b, Igdcc4, Cpxm1, Prrx2
## Six2, Mdk, Gm29478, Kcns1, Foxd1, C130021I20Rik, Cacna2d2, Slc4a10, Fzd7, Erc2
## Negative: Dll4, Rasip1, Exoc3l2, Cd300lg, Arhgef15, Icam2, Mmrn2, Tie1, Cmtm8, Unc45b
## Plvap, Tek, Tead4, Ushbp1, Egfl7, Clic5, Sox18, Kdr, Rassf9, Pcdh12
## Dipk2b, Arl15, Nova2, Myct1, Gpihbp1, Spns2, Galnt15, Pecam1, Kank3, Btnl9
## PC_ 3
## Positive: Arhgap24, Chchd10, Scd1, Gfra1, Srpk3, Cpm, Chst3, Nrxn1, Fcrla, Sbk1
## Gm42836, Spib, Cecr2, 2010309G21Rik, Ctnna3, Gm43388, Ly6d, Cd79a, Bfsp2, Tmem108
## Siglecg, Agbl1, Vpreb3, Il9r, B3gnt5, Tnfrsf13c, Ms4a1, H2-Eb2, Pou2af1, Klhl14
## Negative: Dab2, Gm4951, Rab3il1, Wwp1, Slco2b1, Ophn1, Stab1, Rnase4, Cmklr1, Ang
## Trpv4, F830016B08Rik, Hpgd, Prune2, Serpinb8, Vcam1, Abca9, Gm15635, Mrc1, Ctsb
## Cd36, C3ar1, Reps2, 2610203C22Rik, Fcrls, Lgmn, Hfe, Nfxl1, Ms4a7, Msr1
## PC_ 4
## Positive: Dennd3, Erg, Aff3, 2700054A10Rik, Nxn, Gm16046, Cacna1e, Rgs7, Chrm3, Eya1
## Kcnj15, Gm16070, Ifitm10, Pik3c2b, Tnfsf10, Fcrla, Ly6d, Cd79a, Dchs1, Myzap
## Spib, Siglecg, Trpm6, 2010309G21Rik, Bach2, Pde2a, Gm43388, Pecam1, Gm26760, Dennd5b
## Negative: Kcnj12, Itga7, Ctnna3, Axl, Il34, Sncg, Tagln, Mertk, Nrip2, Rasl12
## Myom1, Synpo2, Akap6, Cap2, Myh11, Myl9, Dgkb, Adap2, Pln, Gm34829
## Jph2, Susd5, Rgs6, Tbx3os1, Casq2, Ppp1r14a, Tpm2, Trarg1, Rgs4, Acta2
## PC_ 5
## Positive: Mmp8, Lcn2, Cd177, Anxa1, Trem3, Dhrs9, B430306N03Rik, 4930438A08Rik, Scrg1, B230208H11Rik
## Cebpe, Acod1, Il1rn, Wfdc21, Camp, Slfn4, Fpr2, F5, Arg2, Fpr1
## Ngp, Clec4d, A530064D06Rik, Trem1, Padi4, Ly6g, Itgb2l, 9830107B12Rik, Cxcr2, Ltf
## Negative: Ptprcap, Gimap7, Ikzf3, Cd2, Gimap3, Fam169b, Gpr174, Bach2, Unc5cl, Ablim1
## Aff3, Slamf6, Sidt1, Myo3b, Tox, Rnf43, Lck, Ctla4, A430093F15Rik, Blk
## Clnk, Btla, Trbc2, Gm39323, Il2rb, Hdac7, Themis, Gprin3, Ctsw, Zap70
# Step 6: Visualize PCA
pca_colors <- c("firebrick1","cyan","gold","blue","black","forestgreen",
"darkorchid1","green","gray","deeppink","chocolate4",
"steelblue","pink","orange")
pca <- DimPlot(pseudo,
reduction = "pca",
group.by = "annotated_clusters",
cols = pca_colors,
pt.size = 3)
pca
pdf(paste0(out, "clustering_QC/pseudobulk_pca.pdf"), width = 6, height = 4)
pca
dev.off()
## png
## 2
# create new object
shiny.obj <- mouse.annotated
VariableFeatures(shiny.obj) <- shiny.obj@assays$SCT@var.features
# set default params
DefaultAssay(shiny.obj) <- "RNA"
Idents(shiny.obj) <- "annotated_clusters"
# create config
names <- colnames(shiny.obj@meta.data)
names <- names[c(22,23,2:21)]
sc.config <- createConfig(obj = shiny.obj,
meta.to.include = names)
# change wd
setwd(out)
# output shiny app folder
makeShinyApp(obj = shiny.obj,
scConf = sc.config,
gene.mapping = TRUE,
shiny.title = "All Clusters")
# manual config edits
sc1conf <- readRDS("shinyApp/sc1conf.rds")
sc1conf[2,4] <- "#B5B9BA|#3385BB|#40BBFF|#a5d5a9|#5dbd64|#1c7e24|#F57C7C|#E42622|#FBB268|#FE8D19|#DE9E83|#A6CEE3|#9D7BBA|#977899"
saveRDS(sc1conf, "shinyApp/sc1conf.rds")